Metrics for Optimizing Searches for Tidally Decaying Exoplanets

To explore the astrophysical properties that support precise transit time estimates, we start with a simplified model for the central time t_c of a transit or eclipse (Carter et al. 2008). This model involves (among other simplifications) neglecting orbital eccentricity and limb-darkening and assuming that the transiting planet is small compared to the star and that the out-of-transit baseline is very accurately estimated. (Numerical experimentation using fully accurate transit light curves shows that this simplified estimate is good to about 10%.) Carter et al. (2008) defined several useful parameters related to the transit:

$$\begin{eqnarray}&&{\tau }_{0}=\displaystyle \frac{{R}_{\star }}{a}\displaystyle \frac{P}{2\pi },\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&b=\displaystyle \frac{a}{{R}_{\star }}\cos i,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&T=2{\tau }_{0}\sqrt{1-{b}^{2}},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\tau =2{\tau }_{0}\left(\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{\star }}\right){\left(1-{b}^{2}\right)}^{-1/2},\end{eqnarray} \tag{ 6 }$$

where R_⋆ is the stellar radius, a is the orbital semimajor axis, P is the orbital period, i is the orbital inclination, b is the impact parameter, T is the total transit duration (defined as the time for the planet's center to cross from limb to limb), and τ is the ingress/egress duration. We also need a transit or eclipse depth δ,

$$\begin{eqnarray}\delta \approx \left\{\begin{array}{ll}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{\star }}\right)}^{2} & \mathrm{transit}\\ {\left(\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{\star }}\right)}^{2}\displaystyle \frac{{I}_{{\rm{p}}}}{{I}_{\star }} & \mathrm{occultation},\end{array}\right.\end{eqnarray} \tag{ 7 }$$

where I_p/⋆ is the planetary/stellar disk-integrated intensity (Winn 2010).

We also need the following parameters as defined in Carter et al. (2008):

$$\begin{eqnarray}&&Q=\sqrt{{\rm{\Gamma }}T}\displaystyle \frac{\delta }{\sigma },\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\theta =\displaystyle \frac{\tau }{T},\end{eqnarray} \tag{ 9 }$$

where Γ is the sampling rate for the transit observations (assumed constant), and σ is the per-point photometric uncertainty. Therefore, Q correlates with total signal-to-noise ratio (S/N) for the transit, and θ is the ratio of the ingress/egress duration to the total transit duration. Based on these definitions, Carter et al. (2008) provided a simplified estimate for the uncertainty on the central time t_c:

$$\begin{eqnarray}&&{\sigma }_{{t}_{{\rm{c}}}}={{TQ}}^{-1}\sqrt{\theta /2}.\end{eqnarray} \tag{ 10 }$$

Plugging in all of the above defined parameters, we find that

$$\begin{eqnarray}&&\begin{array}{l}{\sigma }_{{t}_{{\rm{c}}}}=\sqrt{\tfrac{\tau }{2{\rm{\Gamma }}}}\left(\displaystyle \frac{\sigma }{\delta }\right)=\sigma {{\rm{\Gamma }}}^{-1/2}\\ \ \ \times \ {\left(\displaystyle \frac{{R}_{\star }}{a}\right)}^{1/2}{\left(\displaystyle \frac{P}{2\pi }\right)}^{1/2}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{\star }}\right)}^{-3/2}{\left(1-{b}^{2}\right)}^{-1/4}\\ \ \ \times \ \left\{\begin{array}{ll}1 & \mathrm{transit}\\ {\left(\displaystyle \frac{{I}_{{\rm{p}}}}{{I}_{\star }}\right)}^{-1} & \mathrm{eclipse}\end{array}\right..\end{array}\end{eqnarray} \tag{ 11 }$$

Not surprisingly, the uncertainty increases with the photometric uncertainty and decreases as the transit depth and sampling rate increase. Why, though, does the mid-transit uncertainty increase with ingress/egress duration? Consider a case with a very long ingress/egress (e.g., a near-grazing transit with b → 1), which corresponds to a very nearly V-shaped light curve. In that case, determining the mid-transit time relies on being able to determine when exactly the light curve goes from decreasing with time to increasing with time, with very little transition in between. Without sufficient sampling, for example, the instant of transition would be missed, and the mid-transit time would be highly uncertain.

Figure 1 compares estimates of transit ${\sigma }_{{t}_{{\rm{c}}}}$ for several systems to an estimate for WASP-12 b based on Equation (11). Although it is impossible to estimate the per-point photometric uncertainty σ for any system in general, since the photometric uncertainty depends on the complex details of a particular observation, we can at least include the approximate dependence on stellar magnitude. First, we can relate the photon count rate N to the stellar flux in the bandpass of observation F as N ∝ F. Then we can fold in the relationship between flux and apparent magnitude $m\propto -2.5{\mathrm{log}}_{10}F$. Assuming Poisson statistics gives

$$\begin{eqnarray}&&\sigma \propto {10}^{-m/5}.\end{eqnarray} \tag{ 12 }$$

Zoom In Zoom Out Reset image size

Although simplified, the calculations illustrated in Figure 1 show that there are many systems that we might expect to have smaller timing uncertainties than WASP-12 b and many systems with tidal decay that are expected to be larger. But there are only three with both: WASP-103 b, KELT-16 b, and KELT-1 b. Barros et al. (2022) analyzed combined ground- and space-based transit observations of WASP-103 b, realizing typical timing uncertainties about 50% smaller than those for WASP-12 b reported in Yee et al. (2020). However, Barros et al. (2022) reported no detection of tidal decay but did see tidal deformation of the planet. Likewise, Harre et al. (2023) combined ground- and space-based data for KELT-16 b and found that the BIC favors no tidal decay, but only slightly: BIC = 292.9 for a constant period and BIC = 297.6 for decay. Finally, Baştürk et al. (2023) combined 19 transit observations for the brown dwarf system KELT-1 b and also found no evidence for tidal decay, but they did report possible signs of tidal synchronization of the host star's rotation.

Having developed a sense for the range of timing uncertainties and tidal decay rates, we next turn to how transit observations are transformed into ephemerides, both those that include no decay (linear in the observational epoch E) and those that do include it (quadratic in E).

For a linear fit to a linear ephemeris based solely on transits, we can calculate the uncertainties on T₀ and P using the epoch E for each observed transit and the associated mid-transit time uncertainty σ_t(E). We use σ_t(E) to represent the actually observed uncertainty (as opposed to the analytic uncertainty for a single transit ${\sigma }_{{t}_{{\rm{c}}}}$ or the uncertainty for the predicted future transit time ${\sigma }_{{t}_{\mathrm{tra}}^{\mathrm{pred}}}$). The predicted time and associated uncertainty for the transit time are, respectively,

$$\begin{eqnarray}&&{t}_{\mathrm{tra}}^{\mathrm{pred}}={T}_{0}+{PE},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\sigma }_{{t}_{\mathrm{tra}}^{\mathrm{pred}}}=\sqrt{{\sigma }_{{T}_{0}}^{2}+{E}^{2}{\sigma }_{P}^{2}+2E{\sigma }_{{T}_{0},P}}\approx \sqrt{{\sigma }_{{T}_{0}}^{2}+{E}^{2}{\sigma }_{P}^{2}}.\end{eqnarray} \tag{ 14 }$$

Although the ${\sigma }_{{T}_{0},P}$ term contributes, in practice, it is usually orders of magnitude smaller than the other terms, so we neglect it.

We can estimate the uncertainties analytically using standard linear regression (see Press et al. 2002). First, define

$$\begin{eqnarray}\begin{array}{rcl}S & = & \displaystyle \displaystyle \sum _{E\in \mathrm{transits}}\left(1/{\sigma }_{t(E)}^{2}\right)\\ {S}_{E} & = & \displaystyle \displaystyle \sum _{E\in \mathrm{transits}}\left(E/{\sigma }_{t(E)}^{2}\right)\\ {S}_{{E}^{2}} & = & \displaystyle \displaystyle \sum _{E\in \mathrm{transits}}\left({E}^{2}/{\sigma }_{t(E)}^{2}\right).\end{array}\end{eqnarray} \tag{ 15 }$$

Then

$$\begin{eqnarray}&&{\sigma }_{{T}_{0}}^{2}=\displaystyle \frac{{S}_{{E}^{2}}}{{{SS}}_{{E}^{2}}-{\left({S}_{E}\right)}^{2}}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{\sigma }_{P}^{2}=\displaystyle \frac{S}{{{SS}}_{{E}^{2}}-{\left({S}_{E}\right)}^{2}}.\end{eqnarray} \tag{ 17 }$$

Figure 2 illustrates how adding more and more observations impacts ${\sigma }_{{t}_{\mathrm{tra}}^{\mathrm{pred}}}$ by following the history of transit observations of WASP-12 b.

Zoom In Zoom Out Reset image size

Of course, if we wait for a while before observing another transit, the uncertainty for the next expected transit time will grow as in Equation (14). If we waited long enough, t_wait, that ${\sigma }_{{t}_{\mathrm{tra}}^{\mathrm{pred}}}$ grows beyond some value, then scheduling the next transit observation could be challenging:

$$\begin{eqnarray}&&{t}_{\mathrm{wait}}=\sqrt{{\sigma }_{{t}_{\mathrm{tra}}^{\mathrm{pred}}}^{2}-{\sigma }_{{T}_{0}}^{2}}\left(\displaystyle \frac{P}{{\sigma }_{P}}\right).\end{eqnarray} \tag{ 18 }$$

Systems that have not been observed for t_wait are the ones for which additional transit observations would be most fruitful for improving the linear ephemeris. For the systems considered here, Figure 3 shows the expected time we would have to wait for ${\sigma }_{{t}_{\mathrm{tra}}^{\mathrm{pred}}}$ to grow as large as the transit duration T since the T₀ value reported in the Exoplanet Archive (as of 2023 April 5). Most planets have sufficiently precise linear ephemerides that we would have to wait many years before the uncertainties on their expected ${t}_{\mathrm{tra}}^{\mathrm{pred}}$ grew as large as their transit durations, but uncertainties for a handful are likely large enough to warrant follow-up already, at least based on the Exoplanet Archive data. For example, CoRoT-14 b has an orbital period P = 1.51214 ± 0.00013 days and T₀ = 2,454,787.6702 ± 0.0053 JD (Bonomo et al. 2017), which corresponds to 2008 November. Over the last decade and a half, ${\sigma }_{{t}_{\mathrm{tra}}^{\mathrm{pred}}}$ has grown as large as its transit duration, T = 1.2 hr. Very near the one-to-one line, WASP-103 b was recently observed by the CHEOPS telescope, observations that actually suggest an orbital period increase rather than a decrease (Barros et al. 2022). However, the resulting ephemeris was too recent to have been included in our data, so we do not consider it here. A similar case is TrES-3 b; we did not use more recent observations (e.g., Mannaday et al. 2022) that would likely update its timing uncertainty and increase t_wait. Determining whether individual systems require follow-up or just updated ephemerides is left for future work.

Zoom In Zoom Out Reset image size

Finally, fitting a linear curve to a linear ephemeris would be expected to result in a BIC given by

$$\begin{eqnarray}&&\mathrm{BIC}(\mathrm{lin})=\left(N-2\right)+2\mathrm{ln}N.\end{eqnarray} \tag{ 19 }$$

For a quadratic ephemeris, the predicted time and associated uncertainty for the time of the Eth transit are, respectively,

$$\begin{eqnarray}&&{t}_{\mathrm{tra}}^{\mathrm{pred}}={T}_{0}+{PE}+\displaystyle \frac{1}{2}\left(\displaystyle \frac{{dP}}{{dE}}\right){E}^{2},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{\sigma }_{{t}_{\mathrm{tra}}^{\mathrm{pred}}}\approx \sqrt{{\sigma }_{{T}_{0}}^{2}+{E}^{2}{\sigma }_{P}^{2}+\displaystyle \frac{1}{4}{E}^{4}{\sigma }_{{dP}/{dE}}^{2}},\end{eqnarray} \tag{ 21 }$$

where we have again neglected covariance between fit parameters.

For tidal decay involving a constant phase lag (i.e., a constant value for the star's modified tidal dissipation parameter Q_⋆), dP/dE is given by

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{dP}}{{dE}}=P\displaystyle \frac{{dP}}{{dt}}\approx -\left(26\,\mu {\rm{s}}\ \mathrm{per}\ \mathrm{orbit}\right)\\ \ \ \times \ \left(\displaystyle \frac{{M}_{{\rm{p}}}}{{{\rm{M}}}_{\mathrm{Jup}}}\right){\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{-8/3}{\left(\displaystyle \frac{{R}_{\star }}{{R}_{\odot }}\right)}^{5}\\ \ \ \times \ {\left(\displaystyle \frac{P}{\mathrm{day}}\right)}^{-10/3}{\left(\displaystyle \frac{{Q}_{\star }}{{10}^{5}}\right)}^{-1},\end{array}\end{eqnarray} \tag{ 22 }$$

where M_p is the planetary mass in Jupiter masses (M_Jup = 1.89813 × 10²⁷ kg), M_⋆ is the stellar mass in solar masses (M_⊙ = 1.989 × 10³⁰ kg), R_⋆ is the stellar radius in solar radii (R_⊙ = 6.957 × 10⁸ m), and P is the orbital period in days.

Figure 1 compares estimates of dP/dE for many systems to dP/dE for WASP-12, assuming a WASP-12-like Q_⋆ = 2 × 10⁵. Interestingly, many systems might be expected to exhibit faster tidal decay, and many more systems likely have properties that give rise to more precise transit timing ${\sigma }_{{t}_{{\rm{c}}}}$, at least based on the simplified analytic treatment outlined in Section 2.1.

By analogy with the linear case, we can analytically calculate the uncertainties on the fit parameters. For this calculation, we define

$$\begin{eqnarray}\begin{array}{rcl}{S}_{{E}^{3}} & = & \displaystyle \displaystyle \sum _{E\in \mathrm{transits}}\left({E}^{3}/{\sigma }_{t(E)}^{2}\right),\\ {S}_{{E}^{4}} & = & \displaystyle \displaystyle \sum _{E\in \mathrm{transits}}\left({E}^{4}/{\sigma }_{t(E)}^{2}\right).\end{array}\end{eqnarray} \tag{ 23 }$$

With these definitions,

$$\begin{eqnarray}&&{\sigma }_{{T}_{0}}^{2}=\displaystyle \frac{{S}_{{E}^{3}}^{2}-{S}_{{E}^{2}}{S}_{{E}^{4}}}{{\rm{\Delta }}},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{\sigma }_{P}^{2}=\displaystyle \frac{{S}_{{E}^{2}}^{2}-{{SS}}_{{E}^{4}}}{{\rm{\Delta }}},\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{\sigma }_{{dP}/{dE}}^{2}=\displaystyle \frac{{S}_{E}^{2}-{{SS}}_{{E}^{2}}}{{\rm{\Delta }}},\end{eqnarray} \tag{ 26 }$$

where

$$\begin{eqnarray}&&{\rm{\Delta }}=-{{SS}}_{{E}^{2}}{S}_{{E}^{4}}+{{SS}}_{{E}^{3}}^{2}+{S}_{E}^{2}{S}_{{E}^{4}}-2{S}_{E}{S}_{{E}^{2}}{S}_{{E}^{3}}+{S}_{{E}^{2}}^{3}.\end{eqnarray} \tag{ 27 }$$

Fitting a quadratic curve to a quadratic ephemeris would be expected to result in

$$\begin{eqnarray}&&\mathrm{BIC}(\mathrm{quad})=\left(N-3\right)+3\mathrm{ln}N.\end{eqnarray} \tag{ 28 }$$

We can use these expressions to explore how evidence for tidal decay in the WASP-12 system mounted over the years as a template for finding other systems exhibiting tidal decay. The blue dots in Figure 4 shows the evolution of the tidal decay S/N for WASP-12 b alongside the comparison of the BIC for linear and quadratic fits, ΔBIC ≡ BIC(lin) − BIC(quad). The ΔBIC will grow as the data favor tidal decay. Not surprisingly, as the tidal decay S/N goes up, the BIC preference for the quadratic fit increases, too.

Zoom In Zoom Out Reset image size

Another requirement for the observational constraints on dP/dE to be meaningful is that uncertainties on the linear portion of the ephemeris need to be small compared to the quadratic portion. Otherwise, apparent deviations from a putative linear ephemeris due to tidal decay could be attributed to the uncertainties on the linear ephemeris. This condition translates to

$$\begin{eqnarray}&&{\sigma }_{\mathrm{lin}}\equiv \sqrt{{\sigma }_{{T}_{0}}^{2}+{E}^{2}{\sigma }_{P}^{2}}\lt \displaystyle \frac{1}{2}\left|\displaystyle \frac{{dP}}{{dE}}\right|{E}^{2}.\end{eqnarray} \tag{ 29 }$$

Figure 4 shows how this condition played out for WASP-12 b. The increase in ΔBIC clearly correlates with the growth of the quadratic term in Equation (29).

Considering other systems, Figure 5 shows the cumulative change in orbital period expected due to tidal decay for our collection of systems as compared to the uncertainty on the linear ephemeris. Systems satisfying Equation (29) appear above the orange line. For example, WASP-12, the only system for which tidal decay has been definitively observed, appears above that line, along with several other systems. Several caveats should be considered in evaluating these results, including the fact that we have assumed Q_⋆ = 2 × 10⁵. This is the value inferred for WASP-12, which may exhibit unusually efficient tidal dissipation (Weinberg et al. 2017) and therefore may not be a representative value. For some well-observed systems, such as WASP-18b and WASP-19b, the lack of observed orbital decay to date has been used to constrain their values of Q_⋆ to >10⁷ and >10⁶, respectively (Rosário et al. 2022). Even so, the results point to several systems that merit follow-up transit observations. Several of the systems above the line in Figure 5 have been noted to exhibit period changes. For example, HAT-P-23 has ${\textstyle \tfrac{1}{2}}\left({\textstyle \tfrac{dP}{dE}}\right){E}^{2}=-0.003$ and σ_lin = 0.002 days, while Hagey et al. (2022) reported dP/dt = −5.2 ± 5.8 ms yr⁻¹ (which works out to ΔP ≈ −0.004 days since T₀ for HAT-P-23 b).

Zoom In Zoom Out Reset image size

Finally, we consider the case of fitting a linear curve to a quadratic ephemeris. Analyzing this case is useful because it will allow us to explore how to estimate the BIC thresholds we should look for if we suspect a planet shows signs of tidal decay. (The other combination, fitting a quadratic to a linear ephemeris, would, in principle, result in a quadratic coefficient statistically consistent with zero and the same BIC expression as in Section 2.3.)

To start, consider the linear curve that results from fitting the quadratic ephemeris. Since ∣dP/dE∣ ≪ T₀ and ≪ P (where T₀ and P are the true values for the system), we might suspect that the best-fit values, which we will call ${T}_{0}^{{\prime} }$ and ${P}^{{\prime} }$, would closely resemble the actual values, i.e., ${T}_{0}^{{\prime} }\approx {T}_{0}$ and ${P}^{{\prime} }\approx P$. Indeed, fitting a linear curve to the transit times for WASP-12 b from Yee et al. (2020) returns ${T}_{0}^{{\prime} }$ and ${P}^{{\prime} }$ that match T₀ and P to better than a few parts in 10,000. But the large collection of high-quality data for WASP-12 b means that even this small disagreement is still statistically discrepant. This result comports with the results from Section 2.3; those systems for which we have sufficient data to detect a nonzero dP/dE are also those for which we have very small error bars on T₀ and P. Therefore, in order to calculate the BIC for fitting a linear curve to a quadratic ephemeris, we will need to also calculate ${T}_{0}^{{\prime} }$ and ${P}^{{\prime} }$, which can be written as

$$\begin{eqnarray}&&{T}_{0}^{{\prime} }={T}_{0}+\left(\displaystyle \frac{{dP}}{{dE}}\right){\rm{\Delta }}{T}_{0}^{{\prime} },\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{P}^{{\prime} }=P+\left(\displaystyle \frac{{dP}}{{dE}}\right){\rm{\Delta }}{P}^{{\prime} },\end{eqnarray} \tag{ 31 }$$

where ${\rm{\Delta }}{T}_{0}^{{\prime} }$ and ${\rm{\Delta }}{P}^{{\prime} }$ are the corrections we need to work out. As outlined in the Appendix, standard linear regression gives the following formulae for ${\rm{\Delta }}{T}_{0}^{{\prime} }$ and ${\rm{\Delta }}{P}^{{\prime} }$:

$$\begin{eqnarray}&&{\rm{\Delta }}{T}_{0}^{{\prime} }=\left(\displaystyle \frac{{S}_{{E}^{2}}^{2}-{S}_{{E}^{3}}{S}_{E}}{{S}_{{E}^{2}}S-{S}_{E}^{2}}\right),\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{P}^{{\prime} }=\left(\displaystyle \frac{{S}_{{E}^{3}}S-{S}_{{E}^{2}}{S}_{E}}{{S}_{{E}^{2}}S-{S}_{E}^{2}}\right).\end{eqnarray} \tag{ 33 }$$

Now we can calculate the resulting χ² for fitting a linear curve to a quadratic ephemeris,

$$\begin{eqnarray}\begin{array}{rcl}{\chi }^{2} & = & \displaystyle \displaystyle \sum _{E\in \mathrm{transits}}{\left(\displaystyle \frac{t(E)-{P}^{{\prime} }E-{T}_{0}^{{\prime} }}{{\sigma }_{t(E)}}\right)}^{2}\\ & \approx & \displaystyle \frac{1}{4}{\left(\displaystyle \frac{{dP}}{{dE}}\right)}^{2}\displaystyle \displaystyle \sum _{E\in \mathrm{transits}}{\left(\displaystyle \frac{{E}^{2}-{\rm{\Delta }}{P}^{{\prime} }E-{\rm{\Delta }}{T}_{0}^{{\prime} }}{{\sigma }_{t(E)}}\right)}^{2}+\left(N-2\right),\end{array}\end{eqnarray} \tag{ 34 }$$

and the difference in BIC values for a linear curve and a quadratic curve, both fit to a quadratic ephemeris as determined analytically:

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}\mathrm{BIC} & = & \displaystyle \frac{1}{4}{\left(\displaystyle \frac{{dP}}{{dE}}\right)}^{2}\displaystyle \displaystyle \sum _{E\in \mathrm{transits}}{\left(\displaystyle \frac{{E}^{2}-{\rm{\Delta }}{P}^{{\prime} }E-{\rm{\Delta }}{T}_{0}^{{\prime} }}{{\sigma }_{t(E)}}\right)}^{2}\\ & & -\ \mathrm{ln}N+1.\end{array}\end{eqnarray} \tag{ 35 }$$

Of course, given a set of already observed transits, we could easily calculate the ΔBIC. The benefit of Equation (35) is that we can estimate the ΔBIC expected for a sequence of planned transit observations that have yet to be conducted (given reasonable estimates for the expected timing uncertainties). The dashed orange line in Figure 4 shows how closely the analytic approximation matches the numerical result obtained by directly comparing a linear to a quadratic fit.

First, we consider hypothetical cases to gauge how effectively tidal decay could, in principle, be detected by various observing strategies for planetary systems with definite quadratic ephemerides. For many of the calculations in this section, we assumed a tidal decay rate equal to WASP-12 b's, P/∣dP/dt∣ = 2.3 Myr or dP/dE = 86.7 μs orbit⁻¹ (Yee et al. 2020), and a constant transit timing uncertainty, ${\sigma }_{t(E)}=\mathrm{const}.$ Strictly, ΔBIC depends on the transit epoch E and not on orbital period, but to give a sense of the timescales over which observational campaigns might be conducted, we assumed WASP-12 b's orbital period P = 1.091419649 days to convert from E to years.

To begin with, we consider some overly simple observational campaigns (left panel of Figure 6). (N.B.: Throughout this section, the x- and y-axes of different panels often do not match up.) For the blue and orange lines, that uncertainty was taken as equal to the median for the WASP-12 b data set from Yee et al. (2020), σ_t(E) = 〈σ_W12b〉 = 0.00032 days, while the green line shows the result for an uncertainty that is 10 times larger (0.0032 days). The blue and green lines show how ΔBIC would grow if we could (unrealistically) observe every transit, while the orange line shows what would happen if we observed every 10th transit.

Zoom In Zoom Out Reset image size

As previously stated, ΔBIC > 0 indicates a statistical preference for a quadratic over a linear ephemeris, and all curves in the left panel of Figure 6 show ΔBIC initially dropping from zero into negative values. Intuitively, this behavior reflects the need for curvature in the ephemeris to build up over time so that the quadratic term ($\tfrac{1}{2}| {dP}/{dE}| {E}^{2}$; see Equation (20)) grows sufficiently large that a linear regression is impacted. In other words, we have to wait for a while after a transiting planet is discovered to spot tidal decay. Equation (35) indicates that that crossover point depends on the total number of observations N, the timing/frequency of those observations (the summation term), the timing uncertainty, and the tidal decay rate dP/dE.

The right panel of Figure 6 shows the dependence of the crossover epoch E_crossover for a desired ΔBIC value on the ratio ∣dP/dE∣/σ_t(E), assuming that every transit since a planet's discovery is observed. If, for example, ΔBIC = 10 were the goal of an observing program (dashed orange curve) for a system with WASP-12 b–like properties (∣dP/dE∣/σ_t(E) ≈ 3 × 10⁻⁶), then a minimum of about 1000 transits would need to be observed. Approaching the situation from the opposite direction, a survey including nearly every single transit up to E = 1500 and with WASP-12 b–like timing uncertainties would be expected to achieve ΔBIC ≈ 100 if the system actually exhibited WASP-12 b's tidal decay. In this way, we can apply Equation (35) to a particular observing program to determine a reasonable threshold for decay detection.

Next, we consider somewhat more realistic observing campaigns (Figure 7). The top left panel compares ΔBIC growth for one transit observation and two observations in 1 Earth yr. In about 9 Earth yr, the curve corresponding to twice-annual observations (dashed orange) has grown to nearly twice the ΔBIC for once-annual observations (solid blue). The top right panel compares ΔBIC growth for one observation every 2 months all Earth year-round and the other involving one observation every 2 months for 6 months (dashed orange). Here the two curves weave over one another, suggesting little advantage of one program over the other. This result is not surprising, since little curvature develops in the ephemeris over 6 months.

Zoom In Zoom Out Reset image size

Next, consider the bottom left panel of Figure 7. The solid blue curve involves two consecutive transit observations each Earth year, while the dashed orange curve involves 10 consecutive observations each Earth year. Not surprisingly, ΔBIC increases much more rapidly for the latter program than the former because, in each observing session, significantly more transits are collected. The dashed–dotted green curve involves the same total number of transits as the dashed orange curve but randomly phased (i.e., not necessarily consecutive transits), illustrating that the timing of observations over a short (compared to the decay time) timescale has minimal impact on the evolution. For this particular instantiation, the final ΔBIC ends up slightly below the dashed orange curve, but other examples (not shown) have ΔBIC equal to or even slightly above the dashed orange curve, depending on exactly how the observations are timed.

Finally, consider the bottom right panel of Figure 7. This panel shows how ground-based observations can combine with observations from a TESS-like mission to detect tidal decay. For this calculation, we first assumed that every transit was observed for a WASP-12–like system during a 27 day period in each year, followed by no observations for 25 TESS sectors, and then another sequence of transits were observed, etc. The solid blue line shows this scenario. The dashed orange line shows the same program except with a single ground-based transit included every 6 months. The dashed–dotted green line shows the same program except with six ground-based transit observations randomly spread out during a 6 month observing season. Not surprisingly, the ΔBIC for the TESS + ground programs grows more quickly, demonstrating the power of combining the two approaches. The dashed–dotted green line (six ground-based transits) significantly exceeds the solid blue line once the signal of tidal decay starts to emerge and ΔBIC grows large, while the dotted orange line (one ground-based transit) only modestly exceeds it.

Although the approach here needs to be tailored to each specific observing program for detailed predictions, these results illustrate its general utility. They show that detecting tidal decay requires collecting regular observations and allowing sufficient time for decay to manifest. In general, a significant increase in ΔBIC requires a significant fractional increase in the number of observations. Adding just a few more observations to an already full observing program does not make much difference unless they are judiciously timed. Two encouraging conclusions of these results: (1) an observing program that can only observe a few times a year can still have an impact, and (2) it may be more worthwhile to double the number of candidates, focusing on the planets most likely to exhibit decay, than to double the number of transits observed for a given planet if a program already involves several observations in a year.

Finally, we consider real systems (Figure 8). These examples all involve systems for which possible tidal decay has been reported. We take observations for TrES-1 b, TrES-2 b, and HAT-P-19 b from Hagey et al. (2022) and KELT-9 b from Harre et al. (2023). For the calculations in this section, we take the correct (and variable) transit timing uncertainties and the corresponding orbital periods to convert epoch E to (Earth) years. To extrapolate the ΔBIC evolution forward in time, we assume that observations continue with the same median frequency as before. For example, TrES-2 b has been observed every four orbits, so we assume that same observing cadence going on past 2020.

Zoom In Zoom Out Reset image size

Hagey et al. (2022) analyzed transit times reported in the Exoplanet Transit Database, ³ and, for TrES-1 b, a tidal decay rate dP/dt = −10.9 ± 2.1 ms yr⁻¹ was favored over a constant period by ΔBIC = 9.7. The left panel of Figure 8 shows a good match between the numerical and analytic estimates for ΔBIC, and, assuming the nominal tidal decay rate, the analytic estimate suggests that ΔBIC ought to exceed 50 within the next few years. However, it may take until about 2030 to reach the same level as reported for WASP-12 b in Yee et al. (2020); not surprising, given that WASP-12 b's dP/dt is about twice as large.

Moving next to TrES-2 b, Hagey et al. (2022) estimated dP/dt = −12.6 ± 2.4 ms yr⁻¹ with tidal decay favored at ΔBIC = 8.3. Again, Figure 8 shows a good match between the numerical and analytic estimates (albeit with considerable scatter in the "Numerical" estimate). Again, the smaller dP/dt than WASP-12 b's means that ΔBIC grows more slowly and may not exceed 50 until 2025. Like the TrES-1 data, the TrES-2 observational data show significant statistical fluctuations in ΔBIC; between $E-{E}_{\min }=500$ and 1000, ΔBIC climbed rapidly before settling back toward zero. The dashed orange line shows that such a rapid increase would not have been expected so soon after the planet's discovery. It remains to be seen whether the recent upward tick in ΔBIC seen in the most recent data represents the beginning of a true increase or whether it too is another statistical fluctuation, although the upward tick is consistent with expectations.

For HAT-P-19 b, Hagey et al. (2022) reported dP/dt = −55.2 ± 7.2 ms yr⁻¹ and dP/dE = 606 μs orbit⁻¹, with tidal decay favored at ΔBIC = 8.3. This dP/dE value is almost seven times that for WASP-12 b, and the analytic ΔBIC reflects this, with a value predicted to exceed that for WASP-12 b by 2024. However, the "Numerical" estimate shows considerable nonmonotonicity, reversing direction and sign several times during the observational baseline. The "Numerical" ΔBIC appears to significantly underperform the analytic estimate. In this context, the right panel of Figure 6 is particularly useful. The average timing uncertainty for the HAT-P-19 b data is 0.0007296 days, about twice the average for WASP-12 b (0.00032 days). If all orbits since discovery had been observed, by $E-{E}_{\min }=1000$, we would expect ΔBIC to exceed 100 (dashed–dotted green line in the right panel of Figure 6). Of course, not every orbit of HAT-P-19 b has been observed, but the fact that the "Numerical" ΔBIC does not yet exceed 10 suggests that perhaps the tidal decay reported for HAT-P-19 b is spurious.

Finally, Harre et al. (2023) combined ground-based, Spitzer, TESS, and CHEOPS transits and eclipses of the ultrahot Jupiter KELT-9 b, which orbits a star at the A/B stellar type boundary, and reported a possible decay rate of dP/dt = −24.42 ± 10.66 ms yr⁻¹ with ΔBIC = 8.4 (although the data show a preference, ΔBIC = 13.2, for apsidal precession). Figure 8 shows that the "Numerical" ΔBIC only became positive with the most recent observations before nosing back to zero with the very last observation. The analytic curve suggests that ΔBIC would not have been expected to cross zero until recently anyway and that it might surpass 20 in 2023. Ivshina & Winn (2022) also considered TESS observations of KELT-9 b and found no evidence for decay. Continued monitoring, especially observations of planetary eclipses, seems likely to resolve whether the system actually experiences tidal decay, which would be especially surprising, since A/B stars are not expected to exhibit significant tidal dissipation (Ogilvie 2014).

What to make of all of these comparisons? One key conclusion is that statistical fluctuations in ΔBIC often appear and may falsely hint at tidal decay. Continued, sustained growth in ΔBIC is probably required to confidently report detection of tidal decay. A calculation like that depicted in the right panel of Figure 6 tailored for a specific campaign provides a way of assessing the threshold ΔBIC beyond which tidal decay may be plausible.